#include <iostream>
#include <iomanip>
#include "CTrackMatrix.h"
#include "CFieldExtern.h"
#include "CTrackMp.h"
#include "CTankDtl.h"
#include "CBeam.h"
#include "CSpaceCharge.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/Vector.h"
#include <fstream>
#include <sstream>
#include <vector>
#include <public.h>
using namespace std;
extern int testnlopt();
int main(int argc,char* argv[])           
{
	cout<<"begin"<<endl;
	ifstream infile("test.input");
	if(!infile)
	{
		cout<<"can't open the input file test.input"<<endl;
		exit(0);
	}
	// 临时的，计算用参数
	// ============================
	stringstream ss;
	string       sRow,sWord;
	int          iTmp;
	string       sTmp;
	//=============================
	//整体参数
	//=============================
	int       iLatticeType,iLatticeStart;
	double    dFreq,dEnergyIn,dEnergyOut=0,dPowerLimit=0;
	int       iTankNum,iCellNum,iCellNumLimit=0;
	int       flag1,flag2;
	//=============================
	//space charge参数
	int    nx,ny,nz;
	double dx,dy,dz;
	//=============================
	//DTL cell参数
	//=============================
	double   dLeft,dRight,dSynPhi,dLQ1,dLQ2,dSQ1,dSQ2;
	//=============================
	double   dBasicGradient,dBasicQuadLength,dQuadLength;
	double   dStep=0.02;  //integrition step length
	vector<double>             vdEfield,vdSynPhi;
	vector<vector<double>* >  *vvdGeom;
	CTankDtl                  *DtlTank;
	CTrackMatrix *matrix;
	CTrackMp     *multiParticle;
	vdEfield.resize(0);
	vdSynPhi.resize(0);
	multiParticle=new CTrackMp();
	matrix=new CTrackMatrix();
	flag1=0;flag2=0;
	//line 0
	while(getline(infile,sRow))
	{
		double dEnergy;
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		//这里的twiss参数是所有粒子的还是rms的？初步判断是全粒子的！
		ss>>dEnergy;
		matrix->setEnergyIn(dEnergy);
		//multiparticle  ??when the beam need to generate their own distribution
		break;
	}
	//line 1
	while(getline(infile,sRow))
	{
		double dAlfax,dAlfay,dAlfaz,dBetax,dBetay,dBetaz,dEmitx,dEmity,dEmitz;
		double    dCurrent,dDutyFactor;
		int       iParticleNumber;
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		ss>>dAlfax>>dBetax>>dEmitx>>dAlfay>>dBetay>>dEmity>>dAlfaz>>dBetaz>>dEmitz>>dStep;
		ss>>iParticleNumber>>dCurrent>>dDutyFactor;
		//这里的twiss参数是所有粒子的还是rms的？初步判断是全粒子的！
		multiParticle->beam->setBeam(dAlfax,dAlfay,dAlfaz,dBetax,dBetay,dBetaz,dEmitx,dEmity,dEmitz,iParticleNumber,dCurrent,dDutyFactor,dStep);
		multiParticle->SetDt(dStep);
		multiParticle->dt=dStep*1.0e-9;
    matrix->setCurrent(dCurrent);
		//may also be used in matrix class
		break;
	}
	//line 2
	while(getline(infile,sRow))
	{
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		const  char *aTmp;
		ss>>sTmp;
		aTmp=sTmp.c_str();
		cout<<"begin beam dist reading"<<endl;
		if((aTmp[0]>64&&aTmp[0]<=90)||(aTmp[0]>97&&aTmp[0]<=122))
		{
			cout<<"aTmp[0] is "<<aTmp[0]<<endl;
			multiParticle->beam->beamDistribution(sTmp);
			cout<<"调用beamDistribution程序"<<endl;
		}
		else if(aTmp[0]>48&&aTmp[0]<58)
		{
			ss.clear();
			ss.str(sTmp);
			ss>>iTmp;
			multiParticle->beam->beamDistribution(iTmp);
		}
		else
			cout<<"Distribution 后面的parameter错误"<<endl;
		cout<<"complete Distribution reading"<<endl;
		break;
	}
	//line 3
	while(getline(infile,sRow))
	{
		double dStepMp;
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		else 
		{ 
			ss>>dStepMp; 
			if(dStepMp>0)
			{
				multiParticle->SetDt(dStepMp);
				multiParticle->dt=dStepMp*1.0e-9;
				flag1=1;
			}
			else
				flag1=0;
			break;
		}
	}
	//line 4
	while(getline(infile,sRow))
	{
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		ss>>dFreq;
		multiParticle->SetFreq(dFreq);
		matrix->setFreq(dFreq);
		break;
	}
	//line 5
	while(getline(infile,sRow))
	{
		double dtmp;
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		ss>>dtmp; 
		multiParticle->SetBoreRadius(dtmp);
		break;
	}

	//line 6
	while(getline(infile,sRow))
	{
		double dtmp;
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		ss>>nx>>ny>>nz>>dx>>dy>>dz;
		multiParticle->SC->InGrid(nx,ny,nz,dx,dy,dz);
		multiParticle->SC->formGrid();
		multiParticle->SC->beam=multiParticle->beam;
		multiParticle->SC->setFreq(dFreq);
		break;
	}
	//section 2
	//	getline(infile,sRow);
	multiParticle->getBeam();
	vector<double> vdMI_Twiss(9);  //for matrix input
	vector<double> vdMI_Beam(6);  //for matrix input
	vdMI_Twiss=multiParticle->getMatrixInput_Twiss();
	vdMI_Beam=multiParticle->getMatrixInput_Beam();
	matrix->BeamVector(vdMI_Beam);
	matrix->TwissVector(vdMI_Twiss);

	while(getline(infile,sRow))
	{
		double dtmp1,dtmp2;
		if(sRow==""||sRow==";") continue;
		ss.clear();    ss.str(sRow); ss>>sWord;                 
		if(sWord==""||sWord==";") continue;
		if(sWord=="endFile") break;
		if(sWord=="Drift") 
		{
			ss>>dtmp1>>dtmp2;
			if(flag1==1)
				multiParticle->Drift_Matrix(dtmp1);
			matrix->Drift(dtmp1);
		}
		if(sWord=="RFGap") 
		{
      double dtmp3;
			ss>>dtmp1>>dtmp2>>dtmp3;
			matrix->RFGap(dtmp1,dtmp2,dtmp3);
		}
		if(sWord=="Matrix_Beam") 
		{
			for(int i=0;i<vdMI_Twiss.size();i++)
			{ ss>>dtmp1;
				vdMI_Twiss[i]=dtmp1;
			}
			matrix->TwissVector(vdMI_Twiss);
		}
		if(sWord=="DtlTank_file") 
		{
		  ss>>sTmp;
			matrix->DtlTank_File(sTmp);
			//matrix->DtlTank(sTmp);
		}
		else if(sWord=="Quad")
		{ 
			ss>>dtmp1>>dtmp2;  //L S
			if(flag1==1)
				multiParticle->Quad_Matrix(dtmp1,dtmp2);
			matrix->Quad(dtmp1,dtmp2);
		}
		else if(sWord=="DtlTank")  
		{
			double dETmp,dSynPhiTmp;
			flag2=1;
			ss>>iTankNum>>dFreq>>dETmp>>dSynPhiTmp>>dEnergyIn;
			ss>>dEnergyOut>>dPowerLimit>>iCellNumLimit;       //attention! what if there is only 5 parameters  ?
			vdEfield.push_back(dETmp);
			vdSynPhi.push_back(dSynPhiTmp);
			//Track??????????????!!!!!!!!!!!!!!!!
			DtlTank=new CTankDtl(iTankNum,dFreq,vdEfield,vdSynPhi,dEnergyIn,dEnergyOut,dPowerLimit,iCellNumLimit);
			DtlTank->ReadSfTab("sf");
			DtlTank->Update();
			multiParticle->SetFreq(dFreq);
			vvdGeom=DtlTank->vvdStruct;
			//注意删除时，先删除Dltank,然后matrix，不需要删除vvdGeom,因为没有new,删除DtlTank时候，
			//已经删除了vvdStruct，二者指向同一个地址
		}
		else if(sWord=="QuadFile")  
		{
			ss>>sTmp;
			if(flag1==1)
				multiParticle->QuadForFieldMap(sTmp.c_str());
		}
		else if(sWord=="DTL_FieldMap")  
		{
			double dtmp;
			string sTmp2;
			ss>>sTmp>>sTmp2>>dtmp;
			if(flag1==1)
				multiParticle->DTL_FieldMap(sTmp,sTmp2,dtmp);
		}
	}
	 multiParticle->OutBeamRmsSize();
	multiParticle->beam->DisTest();
	//multiParticle->testParticle();
	multiParticle->OutTwiss();
	//matrix->Drift(1.70712);
	//matrix->getTwiss();
	//matrix->getPhaseAdvance();
	//matrix->getSigMatrix();
//	matrix->Match_MMF();
	//matrix->test();
  //matrix->InitLattice();
	//matrix->Equip_Tradtion();
	//matrix->Equip_Match();
	//matrix->getSigMatrix();
	//testnlopt();

	/*
		 else if(sWord=="Efield")
		 {
		 flag2=1;
		 vdEfield.resize(0);
		 while(ss>>dTmp)
		 vdEfield.push_back(dTmp);
		 }
		 else if(sWord=="SynPhi")
		 {
		 flag3=1;
		 vdSynPhi.resize(0);
		 while(ss>>dTmp)
		 vdSynPhi.push_back(dTmp);
		 }

		 else if(sWord=="Quad")
		 {
		 flag6=1;
		 ss>>iLatticeType>>iLatticeStart>>dBasicGradient>>dBasicQuadLength;
		 }
		 else if(sWord=="Matrix")
		 {
		 flag7=1;
		 matrix=new CTrackMatrix(dFreq,dEnergyIn,dStep);
		 }
		 else if(sWord=="Out1")        //output the geometry information from class CTankDtl
		 {;}
		 else if(sWord=="begin")
		 {
		 multiParticle->beam->setEnergyS(3.0258);
		 double dxm,dxpm,dym,dypm,dzm,dzpm;
	//都是实验室坐标系长度单位是m，dxpm，dypm，dzpm单位m/s
	dxm=(*multiParticle->beam->getRmsSize())[0];
	dym=(*multiParticle->beam->getRmsSize())[1];
	dzm=(*multiParticle->beam->getRmsSize())[2];
	dxpm=(*multiParticle->beam->getRmsSizep())[0];
	dypm=(*multiParticle->beam->getRmsSizep())[1];
	dzpm=(*multiParticle->beam->getRmsSizep())[2];
	CLHEP::HepVector hepvec;
	CLHEP::HepMatrix hepmat;
	hepvec=matrix->BeamVector(dxm,dxpm,dym,dypm,dzm,dzpm);
	hepmat=matrix->Drift(4.1641,3.0258);
	 */
	//     hepvec=(hepvec)*(hepmat);
	//=================multiparticle=================================  
	//  multiParticle->beam->setStep(dStep/3.0);
	//multiParticle->beam->setStep(10e-9);  //the unit of time is second,not nani-second
	// multiParticle->Drift(0);
	/*
		 if(flag8==1)
		 {
		 if(flag1==1)
		 {
		 DtlTank=new CTankDtl(iTankNum,dFreq,vdEfield,vdSynPhi,dEnergyIn,dEnergyOut,dPowerLimit,iCellNumLimit);
		 DtlTank->ReadSfTab("sf");
		 DtlTank->Update();
		 multiParticle->SetFreq(dFreq);
		 vvdGeom=DtlTank->vvdStruct;
		 iCellNum=vvdGeom->size();
		 cout<<"CTankDtl"<<endl;
		 }
		 for(int i=0;i<1;i++)
		 {
	 */
	/*
	//在parmila精度范围内，xrms，yrms完全吻合
	multiParticle->OutCordFile(11);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,50);
	multiParticle->OutCordFile(21);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,-50);
	multiParticle->OutCordFile(31);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,50);
	multiParticle->OutCordFile(41);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,-50);
	multiParticle->OutCordFile(51);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,50);
	multiParticle->OutCordFile(61);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,-50);
	multiParticle->OutCordFile(71);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,50);
	multiParticle->OutCordFile(81);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,-50);
	multiParticle->OutCordFile(91);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,50);
	multiParticle->OutCordFile(101);
	multiParticle->SetDt(0.05);  //the unit of time is second,not nani-second
	multiParticle->Quad(4,-50);
	multiParticle->OutCordFile(111);
	 */
	/*
		 multiParticle->beam->setStep(0.8327e-9);  //the unit of time is second,not nani-second
		 multiParticle->Quad(50,1);
	 */
	//  SC->scheff(dStep/3.0);  //参数是dt
	//  multiParticle->ThinGap(1,90,324,15);
	//ThinGap(double _dE0TL, double _phi, double _dFreqThinGap, double _aperture)

	//       //计算长drift的twiss参数
	//         multiParticle->SetDt(0.5);  //the unit of time is second,not nani-second
	//         multiParticle->OutCordFile(00);
	//         multiParticle->Drift(12.00907);
	//         multiParticle->OutCordFile(10);
	//        // SC->InputBeam(multiParticle->beam,dCurrent);
	//         //SC->scheff(5.0*2.0e-9);  //参数是dt
	//         multiParticle->SetDt(0.5);  //the unit of time is second,not nani-second
	//         multiParticle->Drift(12.00907);
	//         multiParticle->OutCordFile(20);
	//         multiParticle->SetDt(0.5);  //the unit of time is second,not nani-second
	//         multiParticle->Drift(12.00907);
	//         multiParticle->OutCordFile(30);
	//         multiParticle->SetDt(0.5);  //the unit of time is second,not nani-second
	//         multiParticle->Drift(12.00907);
	//         multiParticle->OutCordFile(40);
	//       // SC->scheff(5.0*2.0e-9);  //参数是dt
	//         multiParticle->Drift(0.1200907);
	//         multiParticle->OutCordFile(200);
	//         multiParticle->Drift(0.1200907);
	//         multiParticle->Drift(0.1200907);
	//       // multiParticle->OutCordFile(400);
	//         multiParticle->Drift(0.1200907);
	//       // multiParticle->OutCordFile(500);
	//         multiParticle->Drift(0.1200907);
	//       // multiParticle->OutCordFile(600);
	//        multiParticle->Drift(0.1200907);
	//       // multiParticle->OutCordFile(700);
	//         multiParticle->Drift(0.1200907);
	//       // multiParticle->OutCordFile(800);
	//
	//		    CDtlCellDesign *celldesign=DtlTank->cDtlCell;             //problems just appeared today
	//		    multiParticle->SetCellDesign(celldesign);
	// multiParticle->ThinGap(3.0258,2,90,324,15);
	//without acceleration,2MV/m
	//==========-90deg
	//    multiParticle->beam->setStep(0.83267e-9);  //the unit of time is second,not nani-second
	//    multiParticle->Quad(2,0.5);
	//    multiParticle->beam->setStep(0.83267e-9);  //the unit of time is second,not nani-second
	//    multiParticle->Drift(0.120091);
	/*
	//quad only
	multiParticle->SetPhaseMove(-90);
	multiParticle->OutCordFile(100);
	multiParticle->SetTTF(0.7724,0.0649,0.5045,0.0563);
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2,  5, 2, -5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2, -5, 2,  5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2,  5, 2, -5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2, -5, 2,  5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2,  5, 2, -5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2, -5, 2,  5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2,  5, 2, -5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2, -5, 2,  5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2,  5, 2, -5, 3.0258, 3.0258, 0);                            
	multiParticle->DtlCell(3.70655, 3.70655, 2.2779,  0, -90, 2, -5, 2,  5, 3.0258, 3.0258, 0);                            
	 */
	/*
		 multiParticle->SetPhaseMove(-90);
		 multiParticle->OutCordFile(100);
		 multiParticle->SetTTF(0.7661,0.0649,0.5046,0.0562);
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->OutCordFile(200);
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->OutCordFile(300);
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
		 multiParticle->DtlCell(3.69485, 3.69485, 2.2695,  4, -90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
	 */
	/* 	//multiParticle->SetPhaseMove(0.00);
			multiParticle->SetPhaseMove(-215.00);
			multiParticle->OutCordFile(100);
	//multiParticle->SetTTF(0.7785,0.0648,0.5044,0.0563);
	multiParticle->FieldMap("OUTSF7.TXT");
	//multiParticle->DtlCell(3.7183, 3.7183, 2.2865,  4,  90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
	multiParticle->QuadForFieldMap("Quadinfo.DST");
	multiParticle->Element_FieldMap(40);
	//FileChangeForImpactT_field2d("OUTSF7.TXT.2");
	//multiParticle->DTL_Fourier("line2.txt",1);
	//multiParticle->DTL_Fourier("line3.TXT",1.0);
	//FileChangeForImpactT_partcl("part_rfq.TXT.0mAmactch");
	//multiParticle->DTL_Fourier("line.TXT",0.07);

	//				  multiParticle->Quad_Matrix(2,-50);
	//					multiParticle->Drift_Matrix(6);
	//				  multiParticle->Quad_Matrix(4,50);
	//					multiParticle->Drift_Matrix(6);
	//				  multiParticle->Quad_Matrix(4,-50);
	//					multiParticle->Drift_Matrix(6);
	//				  multiParticle->Quad_Matrix(4,50);
	//					multiParticle->Drift_Matrix(6);
	//				  multiParticle->Quad_Matrix(4,-50);
	//					multiParticle->Drift_Matrix(6);
	//				  multiParticle->Quad_Matrix(4,50);
	//					multiParticle->Drift_Matrix(6);
	//				  multiParticle->Quad_Matrix(4,-50);
	//					multiParticle->Drift_Matrix(6);
	// FileChangeForTrack("equalxy2.TXT");
	//multiParticle->Element_FieldMap(7.4837);
	//					SCTest();
	//		    multiParticle->Element_FieldMap(7.4837);
	//multiParticle->Element_FieldMap(35.0);
	// multiParticle->Element_FieldMap(0.74837);
	// 			 multiParticle->Element_FieldMap(10.0);
	//					 multiParticle->Element_FieldMap(10);
	//multiParticle->Element_FieldMap(850.7298);
	// multiParticle->Element_FieldMap(60);
	// multiParticle->MapForTW("PMI1");
	//   multiParticle->SearchInitPhase();
	//   multiParticle->OutCordFile(200);
	//   multiParticle->DtlCell(3.7183, 3.7183, 2.2865,  4, 90, 2,  0.00, 2, -0.0, 3.0258, 3.0258, 0);                            
	//   multiParticle->OutCordFile(300);
	//==============0deg
	//    multiParticle->SetPhaseMove(0);
	//    multiParticle->DtlCell(3.71175,3.76975,2.5151,2,0,2,0.5,2,-0.5
	//     ,3.0258,3.139,0);
	//    multiParticle->DtlCell(3.7793,3.8382,2.5151,2,0,2,0.5,2,-0.5
	//     ,3.0258,3.139,0);
	//=============-30deg
	 */
	/*
		 multiParticle->SetPhaseMove(-30);
		 multiParticle->OutCordFile(100);
		 multiParticle->SetTTF(0.7706,0.0646,0.5040,0.0564);
		 multiParticle->DtlCell(3.70960,3.81270,2.3176,4.0,-30.0,2.0,50.0,2.0,-50.0 ,3.0258,3.22660,0);
		 multiParticle->OutCordFile(200);
		 multiParticle->SetTTF(0.7734,0.0641,0.5028,0.0567);
		 multiParticle->DtlCell(3.83005,3.93345,2.4057,4.0,-30.0,2.0,-50.0,2.0,-50.0 ,3.2266,3.4346,0);
	 */
	//			multiParticle->OutBeamRmsSize();
	//			multiParticle->OutTwiss();
	//			//  emittance("parmila.txt");
	//			// multiParticle->OutCordFile(1);
	//			multiParticle->beam->calRms();
	//		}
	//	}
	//	multiParticle->beam->calRms();
	//of.close();
	//		}
	if(flag2==1)
		delete DtlTank;
	delete multiParticle;
	delete matrix;
}                               

